





import arcpy
from arcpy.sa import *
from arcpy import env

# import numpy as np

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "GIS_TISSA_alt"
        self.alias = "GIS_TISSA_alt"
        # List of tool classes associated with this toolbox
        self.tools = [GIS_TISSA_csv, GIS_TISSA_GUI]


class GIS_TISSA_csv(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "GIS_TISSA_csv"
        self.description = "GIS_TISSA reading inputs from a .csv file"
        # self.description = "Some description..."
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
		# Thermal camera image in raster format

        param0 = arcpy.Parameter(
		    name="Digital Elevation Model",
            displayName="Digital Elevation Model",
            datatype="DERasterDataset",
            parameterType="Required",
            direction="Input"
        )
        param0.value="dem.tif"
		
        param1 = arcpy.Parameter(
		    name="Soil Classes",
            displayName="Soil Classes",
            datatype="DEFeatureClass",
            parameterType="Required",
            direction="Input"
        )
        param1.value="soils.shp"

        param2 = arcpy.Parameter(
		    name="Vegetation Classes",
            displayName="Vegetation Classes",
            datatype="DEFeatureClass",
            parameterType="Required",
            direction="Input"
        )
        param2.value="trees.shp"

        param3 = arcpy.Parameter(
		    name="Soils CSV file",
            displayName="Soils CSV File",
            datatype="File",
            parameterType="Required",
            direction="Input"
        )
        param3.value="Soil_prop.csv"

        param4 = arcpy.Parameter(
		    name="Vegetation CSV file",
            displayName="Vegetation CSV File",
            datatype="File",
            parameterType="Required",
            direction="Input"
        )
        param4.value="Veg_prop.csv"

        param5 = arcpy.Parameter(
		    name="Minimum slope value",
            displayName="Minimum Slope Value",
            datatype="GPDouble",
            parameterType="Required",
            direction="Input"
        )
        param5.value=5

        param6 = arcpy.Parameter(
		    name="Output geodatabase",
            displayName="Output Geodataset",
            datatype="DEWorkspace",
            parameterType="Required",
            direction="Input"
        )

        param7 = arcpy.Parameter(
		    name="DEM elevation standard deviation",
            displayName="DEM elevation sTandard deviation",
            datatype="GPDouble",
            parameterType="Required",
            direction="Input"
        )
        param7.value=2

        param8 = arcpy.Parameter(
		    name="Expected seismic ground acceleration",
            displayName="Expected seismic ground acceleration",
            datatype="GPDouble",
            parameterType="Required",
            direction="Input"
        )
        param8.value=0.39
		

        param11 = arcpy.Parameter(
		    name="Slope_method",
            displayName="Slope calculation method",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param11.filter.type="ValueList"
        param11.filter.list=["second-order accurate central difference approximation (PISA-m method)", "Eight neighbour pixel method (ArcMap method)"]
        param11.value="Eight neighbour pixel method (ArcMap method)"
		
        params = [param0, param1, param2, param3, param4, param5, param6, param7, param8, param11]
				  
        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def execute(self, parameters, messages):
        """The source code of the tool."""
        arcpy.CheckOutExtension("Spatial")
        arcpy.CheckOutExtension("3D")

        arcpy.AddMessage("Assigning inputs")
		
		# Assigning inputs
        DEM_1 = parameters[0].valueAsText
        Soil_map = parameters[1].valueAsText
        Veg_map = parameters[2].valueAsText
        Soil_CSV = parameters[3].valueAsText
        Veg_CSV = parameters[4].valueAsText
        Minimum_slope_value = parameters[5].valueAsText
        wd=parameters[6].valueAsText
        DEM_SD = parameters[7].valueAsText
        SeisA0 = parameters[8].valueAsText
        SlopeMeth = parameters[9].valueAsText
        
        DEM_2=Raster(DEM_1)
		
		
		# Define workspace and overwritting permission
        arcpy.env.workspace=wd
        arcpy.env.overwriteOutput = True

		# Extract the pixel size and raste extent from the DEM
        cellsz=arcpy.GetRasterProperties_management(DEM_2,"CELLSIZEX")
        x1=arcpy.GetRasterProperties_management(DEM_1,"LEFT")
        x2=arcpy.GetRasterProperties_management(DEM_1,"RIGHT")
        y1=arcpy.GetRasterProperties_management(DEM_1,"BOTTOM")
        y2=arcpy.GetRasterProperties_management(DEM_1,"TOP")

		#Create constant raster with water weight values
        gw=CreateConstantRaster(9810,"Float",cellsz,Extent(x1, y1, x2, y2))
		

        if SlopeMeth=="second-order accurate central difference approximation (PISA-m method)":

            arcpy.AddMessage("Slope calculation compatible with PISA-m")
		
		    # Slope calculation compatible with PISA-m
            cellsz2=cellsz.getOutput(0)
            cellsz4=float(cellsz2)
            nshift="-"+cellsz.getOutput(0)
		
            arcpy.Shift_management(DEM_2,"DEM_N",0,cellsz)
            arcpy.Shift_management(DEM_2,"DEM_S",0,nshift)
            arcpy.Shift_management(DEM_2,"DEM_E",cellsz,0)
            arcpy.Shift_management(DEM_2,"DEM_W",nshift,0)
		
            DEM_N2=Raster("DEM_N")
            DEM_S2=Raster("DEM_S")
            DEM_E2=Raster("DEM_E")
            DEM_W2=Raster("DEM_W")

            beta=ATan(SquareRoot((DEM_N2-DEM_S2)**2+(DEM_E2-DEM_W2)**2)/(2*cellsz4))
            beta_SD=SquareRoot((8*cellsz4**2*float(DEM_SD)**2)/((4*cellsz4**2+(DEM_N2-DEM_S2)**2+(DEM_E2-DEM_W2)**2)**2))
			
            # arcpy.Delete_management("DEM_N")
            # arcpy.Delete_management("DEM_S")
            # arcpy.Delete_management("DEM_E")
            # arcpy.Delete_management("DEM_W")
			
			
        elif SlopeMeth=="Eight neighbour pixel method (ArcMap method)":
		
            arcpy.AddMessage("Standard ArcMap slope calculation")
			
            beta=Slope(DEM_2,"DEGREE")*3.1415926/180		
			
		    # Slope calculation compatible with PISA-m
            cellsz2=cellsz.getOutput(0)
            cellsz4=float(cellsz2)
            nshift="-"+cellsz.getOutput(0)
		
            arcpy.Shift_management(DEM_2,"DEM_z1",cellsz,cellsz)
            arcpy.Shift_management(DEM_2,"DEM_z2",cellsz,0)
            arcpy.Shift_management(DEM_2,"DEM_z3",cellsz,nshift)
            arcpy.Shift_management(DEM_2,"DEM_z4",nshift,cellsz)
            arcpy.Shift_management(DEM_2,"DEM_z5",nshift,0)
            arcpy.Shift_management(DEM_2,"DEM_z6",nshift,nshift)
            arcpy.Shift_management(DEM_2,"DEM_z7",0,nshift)
            arcpy.Shift_management(DEM_2,"DEM_z8",0,cellsz)
		
            DEM_Z1=Raster("DEM_z1")
            DEM_Z2=Raster("DEM_z2")
            DEM_Z3=Raster("DEM_z3")
            DEM_Z4=Raster("DEM_z4")
            DEM_Z5=Raster("DEM_z5")
            DEM_Z6=Raster("DEM_z6")
            DEM_Z7=Raster("DEM_z7")
            DEM_Z8=Raster("DEM_z8")
			
            ADEM=((DEM_Z1+2*DEM_Z2+DEM_Z3)-(DEM_Z4+2*DEM_Z5+DEM_Z6))/(8*cellsz4)
            BDEM=((DEM_Z6+2*DEM_Z7+DEM_Z3)-(DEM_Z4+2*DEM_Z8+DEM_Z1))/(8*cellsz4)
			
            beta_SD=SquareRoot(3*float(DEM_SD)**2/(8*cellsz4**2)*(ADEM**2+BDEM**2)/((1+ADEM**2+BDEM**2)*SquareRoot(ADEM**2+BDEM**2)))
			
			
            # arcpy.Delete_management("DEM_z1")
            # arcpy.Delete_management("DEM_z2")
            # arcpy.Delete_management("DEM_z3")
            # arcpy.Delete_management("DEM_z4")
            # arcpy.Delete_management("DEM_z5")
            # arcpy.Delete_management("DEM_z6")
            # arcpy.Delete_management("DEM_z7")
            # arcpy.Delete_management("DEM_z8")
			
		
        arcpy.AddMessage("Convert soils and vegetation maps to rasters")
		
		# Convert soils and vegetation maps to rasters
        arcpy.PolygonToRaster_conversion(Soil_map,"SoilInt","Soil_2","","",cellsz)
        arcpy.PolygonToRaster_conversion(Veg_map,"VegInt","Veg_2","","",cellsz)
        arcpy.TableToTable_conversion(Soil_CSV,wd,"soils_table")		
        arcpy.TableToTable_conversion(Veg_CSV,wd,"veg_table")		
		
        arcpy.AddMessage("Create zero valued rasters for all variables and corresponding SD's")
		
		# Create zero valued rasters for all variables and corresponding SD's to be populated from csv tables
        phi4=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        Cs=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        gm=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        gs=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        D=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        Hw=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        Cr=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        qt=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
		
        phi4_SD=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        Cs_SD=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        gm_SD=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        gs_SD=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        D_SD=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        Hw_SD=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        Cr_SD=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
        qt_SD=CreateConstantRaster(0,"Float",cellsz,Extent(x1, y1, x2, y2))
		
        SeisA=CreateConstantRaster(SeisA0,"Float",cellsz,Extent(x1, y1, x2, y2))

        arcpy.AddMessage("Calculate phi and its SD")

		# Calculate phi and its SD
        t2=arcpy.SearchCursor(wd+"\\soils_table")
        for row2 in t2:
            if row2.getValue("param")=="phi" and row2.getValue("dist")=="normal":
                phi0=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat1"))
                phi2=phi4+phi0
                phi4=phi2
                phi0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat2"))
                phi2_SD=phi4_SD+phi0_SD
                phi4_SD=phi2_SD
            elif row2.getValue("param")=="phi" and row2.getValue("dist")=="uniform":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2"))/2
                s_1=(row2.getValue("stat2")-row2.getValue("stat1"))**2/12
                phi0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                phi2=phi4+phi0
                phi4=phi2
                phi0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                phi2_SD=phi4_SD+phi0_SD
                phi4_SD=phi2_SD
            elif row2.getValue("param")=="phi" and row2.getValue("dist")=="Gumbel":
                mean_1=row2.getValue("stat1")+0.577216*row2.getValue("stat2")
                s_1=(3.1415926**2)*(row2.getValue("stat2"))**2/6
                phi0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                phi2=phi4+phi0
                phi4=phi2
                phi0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                phi2_SD=phi4_SD+phi0_SD
                phi4_SD=phi2_SD
            elif row2.getValue("param")=="phi" and row2.getValue("dist")=="triangular":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2")+row2.getValue("stat3"))/3
                s_1=(row2.getValue("stat1")*(row2.getValue("stat1")-row2.getValue("stat3"))+row2.getValue("stat3")*(row2.getValue("stat3")-row2.getValue("stat2"))+row2.getValue("stat2")*(row2.getValue("stat2")-row2.getValue("stat1")))/18
                phi0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                phi2=phi4+phi0
                phi4=phi2
                phi0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                phi2_SD=phi4_SD+phi0_SD
                phi4_SD=phi2_SD
            elif row2.getValue("param")=="phi" and row2.getValue("dist")=="beta-pert":
                mean_1=(row2.getValue("stat1")+4*row2.getValue("stat2")+row2.getValue("stat3"))/6
                s_1=((row2.getValue("stat3")-row2.getValue("stat1"))**2)/36
                phi0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                phi2=phi4+phi0
                phi4=phi2
                phi0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                phi2_SD=phi4_SD+phi0_SD
                phi4_SD=phi2_SD

        arcpy.AddMessage("Calculate Cs and its SD")
				
		# Calculate Cs and its SD
        t2=arcpy.SearchCursor(wd+"\\soils_table")
        for row2 in t2:
            if row2.getValue("param")=="c" and row2.getValue("dist")=="normal":
                Cs0=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat1"))
                Cs2=Cs+Cs0
                Cs=Cs2
                Cs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat2"))
                Cs2_SD=Cs_SD+Cs0_SD
                Cs_SD=Cs2_SD
            elif row2.getValue("param")=="c" and row2.getValue("dist")=="uniform":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2"))/2
                s_1=(row2.getValue("stat2")-row2.getValue("stat1"))**2/12
                Cs0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                Cs2=Cs+Cs0
                Cs=Cs2
                Cs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                Cs2_SD=Cs_SD+Cs0_SD
                Cs_SD=Cs2_SD
            elif row2.getValue("param")=="c" and row2.getValue("dist")=="Gumbel":
                mean_1=row2.getValue("stat1")+0.577216*row2.getValue("stat2")
                s_1=(3.1415926**2)*(row2.getValue("stat2"))**2/6
                Cs0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                Cs2=Cs+Cs0
                Cs=Cs2
                Cs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                Cs2_SD=Cs_SD+Cs0_SD
                Cs_SD=Cs2_SD
            elif row2.getValue("param")=="c" and row2.getValue("dist")=="triangular":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2")+row2.getValue("stat3"))/3
                s_1=(row2.getValue("stat1")*(row2.getValue("stat1")-row2.getValue("stat3"))+row2.getValue("stat3")*(row2.getValue("stat3")-row2.getValue("stat2"))+row2.getValue("stat2")*(row2.getValue("stat2")-row2.getValue("stat1")))/18
                Cs0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                Cs2=Cs+Cs0
                Cs=Cs2
                Cs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                Cs2_SD=Cs_SD+Cs0_SD
                Cs_SD=Cs2_SD
            elif row2.getValue("param")=="c" and row2.getValue("dist")=="beta-pert":
                mean_1=(row2.getValue("stat1")+4*row2.getValue("stat2")+row2.getValue("stat3"))/6
                s_1=((row2.getValue("stat3")-row2.getValue("stat1"))**2)/36
                Cs0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                Cs2=Cs+Cs0
                Cs=Cs2
                Cs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                Cs2_SD=Cs_SD+Cs0_SD
                Cs_SD=Cs2_SD

        arcpy.AddMessage("Calculate gm and its SD")
				
		# Calculate gm and its SD
        t2=arcpy.SearchCursor(wd+"\\soils_table")
        for row2 in t2:
            if row2.getValue("param")=="muw" and row2.getValue("dist")=="normal":
                gm0=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat1"))
                gm2=gm+gm0
                gm=gm2
                gm0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat2"))
                gm2_SD=gm_SD+gm0_SD
                gm_SD=gm2_SD
            elif row2.getValue("param")=="muw" and row2.getValue("dist")=="uniform":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2"))/2
                s_1=(row2.getValue("stat2")-row2.getValue("stat1"))**2/12
                gm0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                gm2=gm+gm0
                gm=gm2
                gm0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                gm2_SD=gm_SD+gm0_SD
                gm_SD=gm2_SD
            elif row2.getValue("param")=="muw" and row2.getValue("dist")=="Gumbel":
                mean_1=row2.getValue("stat1")+0.577216*row2.getValue("stat2")
                s_1=(3.1415926**2)*(row2.getValue("stat2"))**2/6
                gm0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                gm2=gm+gm0
                gm=gm2
                gm0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                gm2_SD=gm_SD+gm0_SD
                gm_SD=gm2_SD
            elif row2.getValue("param")=="muw" and row2.getValue("dist")=="triangular":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2")+row2.getValue("stat3"))/3
                s_1=(row2.getValue("stat1")*(row2.getValue("stat1")-row2.getValue("stat3"))+row2.getValue("stat3")*(row2.getValue("stat3")-row2.getValue("stat2"))+row2.getValue("stat2")*(row2.getValue("stat2")-row2.getValue("stat1")))/18
                gm0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                gm2=gm+gm0
                gm=gm2
                gm0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                gm2_SD=gm_SD+gm0_SD
                gm_SD=gm2_SD
            elif row2.getValue("param")=="muw" and row2.getValue("dist")=="beta-pert":
                mean_1=(row2.getValue("stat1")+4*row2.getValue("stat2")+row2.getValue("stat3"))/6
                s_1=((row2.getValue("stat3")-row2.getValue("stat1"))**2)/36
                gm0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                gm2=gm+gm0
                gm=gm2
                gm0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                gm2_SD=gm_SD+gm0_SD
                gm_SD=gm2_SD

        arcpy.AddMessage("Calculate gs and its SD")
				
		# Calculate gs and its SD
        t2=arcpy.SearchCursor(wd+"\\soils_table")
        for row2 in t2:
            if row2.getValue("param")=="suw" and row2.getValue("dist")=="normal":
                gs0=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat1"))
                gs2=gs+gs0
                gs=gs2
                gs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat2"))
                gs2_SD=gs_SD+gs0_SD
                gs_SD=gs2_SD
            elif row2.getValue("param")=="suw" and row2.getValue("dist")=="uniform":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2"))/2
                s_1=(row2.getValue("stat2")-row2.getValue("stat1"))**2/12
                gs0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                gs2=gs+gs0
                gs=gs2
                gs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                gs2_SD=gs_SD+gs0_SD
                gs_SD=gs2_SD
            elif row2.getValue("param")=="suw" and row2.getValue("dist")=="Gumbel":
                mean_1=row2.getValue("stat1")+0.577216*row2.getValue("stat2")
                s_1=(3.1415926**2)*(row2.getValue("stat2"))**2/6
                gs0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                gs2=gs+gs0
                gs=gs2
                gs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                gs2_SD=gs_SD+gs0_SD
                gs_SD=gs2_SD
            elif row2.getValue("param")=="suw" and row2.getValue("dist")=="triangular":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2")+row2.getValue("stat3"))/3
                s_1=(row2.getValue("stat1")*(row2.getValue("stat1")-row2.getValue("stat3"))+row2.getValue("stat3")*(row2.getValue("stat3")-row2.getValue("stat2"))+row2.getValue("stat2")*(row2.getValue("stat2")-row2.getValue("stat1")))/18
                gs0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                gs2=gs+gs0
                gs=gs2
                gs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                gs2_SD=gs_SD+gs0_SD
                gs_SD=gs2_SD
            elif row2.getValue("param")=="suw" and row2.getValue("dist")=="beta-pert":
                mean_1=(row2.getValue("stat1")+4*row2.getValue("stat2")+row2.getValue("stat3"))/6
                s_1=((row2.getValue("stat3")-row2.getValue("stat1"))**2)/36
                gs0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                gs2=gs+gs0
                gs=gs2
                gs0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                gs2_SD=gs_SD+gs0_SD
                gs_SD=gs2_SD

        arcpy.AddMessage("Calculate D and its SD")
				
		# Calculate D and its SD
        t2=arcpy.SearchCursor(wd+"\\soils_table")
        for row2 in t2:
            if row2.getValue("param")=="D" and row2.getValue("dist")=="normal":
                D0=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat1"))
                D2=D+D0
                D=D2
                D0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat2"))
                D2_SD=D_SD+D0_SD
                D_SD=D2_SD
            elif row2.getValue("param")=="D" and row2.getValue("dist")=="uniform":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2"))/2
                s_1=(row2.getValue("stat2")-row2.getValue("stat1"))**2/12
                D0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                D2=D+D0
                D=D2
                D0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                D2_SD=D_SD+D0_SD
                D_SD=D2_SD
            elif row2.getValue("param")=="D" and row2.getValue("dist")=="Gumbel":
                mean_1=row2.getValue("stat1")+0.577216*row2.getValue("stat2")
                s_1=(3.1415926**2)*(row2.getValue("stat2"))**2/6
                D0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                D2=D+D0
                D=D2
                D0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                D2_SD=D_SD+D0_SD
                D_SD=D2_SD
            elif row2.getValue("param")=="D" and row2.getValue("dist")=="triangular":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2")+row2.getValue("stat3"))/3
                s_1=(row2.getValue("stat1")*(row2.getValue("stat1")-row2.getValue("stat3"))+row2.getValue("stat3")*(row2.getValue("stat3")-row2.getValue("stat2"))+row2.getValue("stat2")*(row2.getValue("stat2")-row2.getValue("stat1")))/18
                D0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                D2=D+D0
                D=D2
                D0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                D2_SD=D_SD+D0_SD
                D_SD=D2_SD
            elif row2.getValue("param")=="D" and row2.getValue("dist")=="beta-pert":
                mean_1=(row2.getValue("stat1")+4*row2.getValue("stat2")+row2.getValue("stat3"))/6
                s_1=((row2.getValue("stat3")-row2.getValue("stat1"))**2)/36
                D0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                D2=D+D0
                D=D2
                D0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                D2_SD=D_SD+D0_SD
                D_SD=D2_SD

        arcpy.AddMessage("Calculate Hw and its SD")
				
		# Calculate Hw and its SD
        t2=arcpy.SearchCursor(wd+"\\soils_table")
        for row2 in t2:
            if row2.getValue("param")=="Hw" and row2.getValue("dist")=="normal":
                Hw0=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat1"))
                Hw2=Hw+Hw0
                Hw=Hw2
                Hw0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*float(row2.getValue("stat2"))
                Hw2_SD=Hw_SD+Hw0_SD
                Hw_SD=Hw2_SD
            elif row2.getValue("param")=="Hw" and row2.getValue("dist")=="uniform":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2"))/2
                s_1=(row2.getValue("stat2")-row2.getValue("stat1"))**2/12
                Hw0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                Hw2=Hw+Hw0
                Hw=Hw2
                Hw0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                Hw2_SD=Hw_SD+Hw0_SD
                Hw_SD=Hw2_SD
            elif row2.getValue("param")=="Hw" and row2.getValue("dist")=="Gumbel":
                mean_1=row2.getValue("stat1")+0.577216*row2.getValue("stat2")
                s_1=(3.1415926**2)*(row2.getValue("stat2"))**2/6
                Hw0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                Hw2=Hw+Hw0
                Hw=Hw2
                Hw0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                Hw2_SD=Hw_SD+Hw0_SD
                Hw_SD=Hw2_SD
            elif row2.getValue("param")=="Hw" and row2.getValue("dist")=="triangular":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2")+row2.getValue("stat3"))/3
                s_1=(row2.getValue("stat1")*(row2.getValue("stat1")-row2.getValue("stat3"))+row2.getValue("stat3")*(row2.getValue("stat3")-row2.getValue("stat2"))+row2.getValue("stat2")*(row2.getValue("stat2")-row2.getValue("stat1")))/18
                Hw0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                Hw2=Hw+Hw0
                Hw=Hw2
                Hw0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                Hw2_SD=Hw_SD+Hw0_SD
                Hw_SD=Hw2_SD
            elif row2.getValue("param")=="Hw" and row2.getValue("dist")=="beta-pert":
                mean_1=(row2.getValue("stat1")+4*row2.getValue("stat2")+row2.getValue("stat3"))/6
                s_1=((row2.getValue("stat3")-row2.getValue("stat1"))**2)/36
                Hw0=(Raster("Soil_2")==row2.getValue("Unit"))*mean_1
                Hw2=Hw+Hw0
                Hw=Hw2
                Hw0_SD=(Raster("Soil_2")==row2.getValue("Unit"))*s_1
                Hw2_SD=Hw_SD+Hw0_SD
                Hw_SD=Hw2_SD


				
        arcpy.AddMessage("Calculate Cr and its SD")

		# Calculate Cr and its SD
        t2=arcpy.SearchCursor(wd+"\\veg_table")
        for row2 in t2:
            if row2.getValue("param")=="Cr" and row2.getValue("dist")=="normal":
                Cr0=(Raster("Veg_2")==row2.getValue("Unit"))*float(row2.getValue("stat1"))
                Cr2=Cr+Cr0
                Cr=Cr2
                Cr0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*float(row2.getValue("stat2"))
                Cr2_SD=Cr_SD+Cr0_SD
                Cr_SD=Cr2_SD
            elif row2.getValue("param")=="Cr" and row2.getValue("dist")=="uniform":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2"))/2
                s_1=(row2.getValue("stat2")-row2.getValue("stat1"))**2/12
                Cr0=(Raster("Veg_2")==row2.getValue("Unit"))*mean_1
                Cr2=Cr+Cr0
                Cr=Cr2
                Cr0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*s_1
                Cr2_SD=Cr_SD+Cr0_SD
                Cr_SD=Cr2_SD
            elif row2.getValue("param")=="Cr" and row2.getValue("dist")=="Gumbel":
                mean_1=row2.getValue("stat1")+0.577216*row2.getValue("stat2")
                s_1=(3.1415926**2)*(row2.getValue("stat2"))**2/6
                Cr0=(Raster("Veg_2")==row2.getValue("Unit"))*mean_1
                Cr2=Cr+Cr0
                Cr=Cr2
                Cr0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*s_1
                Cr2_SD=Cr_SD+Cr0_SD
                Cr_SD=Cr2_SD
            elif row2.getValue("param")=="Cr" and row2.getValue("dist")=="triangular":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2")+row2.getValue("stat3"))/3
                s_1=(row2.getValue("stat1")*(row2.getValue("stat1")-row2.getValue("stat3"))+row2.getValue("stat3")*(row2.getValue("stat3")-row2.getValue("stat2"))+row2.getValue("stat2")*(row2.getValue("stat2")-row2.getValue("stat1")))/18
                Cr0=(Raster("Veg_2")==row2.getValue("Unit"))*mean_1
                Cr2=Cr+Cr0
                Cr=Cr2
                Cr0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*s_1
                Cr2_SD=Cr_SD+Cr0_SD
                Cr_SD=Cr2_SD
            elif row2.getValue("param")=="Cr" and row2.getValue("dist")=="beta-pert":
                mean_1=(row2.getValue("stat1")+4*row2.getValue("stat2")+row2.getValue("stat3"))/6
                s_1=((row2.getValue("stat3")-row2.getValue("stat1"))**2)/36
                Cr0=(Raster("Veg_2")==row2.getValue("Unit"))*mean_1
                Cr2=Cr+Cr0
                Cr=Cr2
                Cr0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*s_1
                Cr2_SD=Cr_SD+Cr0_SD
                Cr_SD=Cr2_SD
				
        arcpy.AddMessage("Calculate qt and its SD")
				
		# Calculate qt and its SD
        t2=arcpy.SearchCursor(wd+"\\veg_table")
        for row2 in t2:
            if row2.getValue("param")=="qt" and row2.getValue("dist")=="normal":
                qt0=(Raster("Veg_2")==row2.getValue("Unit"))*float(row2.getValue("stat1"))
                qt2=qt+qt0
                qt=qt2
                qt0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*float(row2.getValue("stat2"))
                qt2_SD=qt_SD+qt0_SD
                qt_SD=qt2_SD
            elif row2.getValue("param")=="qt" and row2.getValue("dist")=="uniform":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2"))/2
                s_1=(row2.getValue("stat2")-row2.getValue("stat1"))**2/12
                qt0=(Raster("Veg_2")==row2.getValue("Unit"))*mean_1
                qt2=qt+qt0
                qt=qt2
                qt0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*s_1
                qt2_SD=qt_SD+qt0_SD
                qt_SD=qt2_SD
            elif row2.getValue("param")=="qt" and row2.getValue("dist")=="Gumbel":
                mean_1=row2.getValue("stat1")+0.577216*row2.getValue("stat2")
                s_1=(3.1415926**2)*(row2.getValue("stat2"))**2/6
                qt0=(Raster("Veg_2")==row2.getValue("Unit"))*mean_1
                qt2=qt+qt0
                qt=qt2
                qt0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*s_1
                qt2_SD=qt_SD+qt0_SD
                qt_SD=qt2_SD
            elif row2.getValue("param")=="qt" and row2.getValue("dist")=="triangular":
                mean_1=(row2.getValue("stat1")+row2.getValue("stat2")+row2.getValue("stat3"))/3
                s_1=(row2.getValue("stat1")*(row2.getValue("stat1")-row2.getValue("stat3"))+row2.getValue("stat3")*(row2.getValue("stat3")-row2.getValue("stat2"))+row2.getValue("stat2")*(row2.getValue("stat2")-row2.getValue("stat1")))/18
                qt0=(Raster("Veg_2")==row2.getValue("Unit"))*mean_1
                qt2=qt+qt0
                qt=qt2
                qt0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*s_1
                qt2_SD=qt_SD+qt0_SD
                qt_SD=qt2_SD
            elif row2.getValue("param")=="qt" and row2.getValue("dist")=="beta-pert":
                mean_1=(row2.getValue("stat1")+4*row2.getValue("stat2")+row2.getValue("stat3"))/6
                s_1=((row2.getValue("stat3")-row2.getValue("stat1"))**2)/36
                qt0=(Raster("Veg_2")==row2.getValue("Unit"))*mean_1
                qt2=qt+qt0
                qt=qt2
                qt0_SD=(Raster("Veg_2")==row2.getValue("Unit"))*s_1
                qt2_SD=qt_SD+qt0_SD
                qt_SD=qt2_SD

				
        phi=phi4*3.1415926/180
        # phi.save("phi")
        phi_SD=phi4_SD*3.1415926/180
        # phi_SD.save("phi_SD")

        arcpy.AddMessage("Calculating Mean Factor of Safety ")
		
        FS=(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))
        FS.save("FoS_Mean") #static mean
        arcpy.AddMessage("Mean Factor of Safety Calculated")

        arcpy.AddMessage("Calculating Standard Deviation of Factor of Safety")
		
        Dcr=1/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))
        Dcs=1/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))
        Dqt=((Cos(beta))**2*Tan(phi)*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-Sin(beta)*Cos(beta)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        Dgm=(((Cos(beta))**2*Tan(phi)*(D-Hw*D))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-Sin(beta)*Cos(beta)*(D-Hw*D)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        Dgs=(((Cos(beta))**2*Tan(phi)*(Hw*D))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-Sin(beta)*Cos(beta)*(Hw*D)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        DD=(((Cos(beta))**2*Tan(phi)*(gm+(gs-gw-gm)*Hw))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-Sin(beta)*Cos(beta)*(gm+(gs-gm)*Hw)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        DHw=(((gs-gw-gm)*D*(Cos(beta))**2*Tan(phi))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-(gs-gm)*D*Sin(beta)*Cos(beta)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        Dbeta=((2*Cos(beta)*(-Sin(beta))*(qt+gm*D+(gs-gw-gm)*Hw*D)*Tan(phi))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-(qt+gm*D+(gs-gm)*Hw*D)*((Cos(beta))**2+(-Sin(beta))*Sin(beta))*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        Dphi=((qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*(1/(Cos(phi)**2)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))
		
        FS_SD=SquareRoot(Dcr**2*Cr_SD**2+Dcs**2*Cs_SD**2+Dqt**2*qt_SD**2+Dgm**2*gm_SD**2+Dgs**2*gs_SD**2+DD**2*D_SD**2+DHw**2*Hw_SD**2+Dbeta**2*beta_SD**2+Dphi**2*phi_SD**2)
        FS_SD.save("FoS_SD")	# static sd	
        arcpy.AddMessage("Standard Deviation of Factor of Safety Created")

        arcpy.AddMessage("Calculating Probability of FoS < 1")
		
		# Following Abramowitz and Stegun (1964), equation 7.1.26 (pp 299), with an upper error bound of 1.5x10^-7:
        x=-1*(Ln(FS)-0.5*Ln((FS_SD/FS)**2+1))/(SquareRoot(2)*SquareRoot(Ln((FS_SD/FS)**2+1)))
        t=1/(1+0.3275911*Abs(x))
        Erf_1=(1-(0.254829592*t-0.284496736*t**2+1.421413741*t**3-1.453152027*t**4+1.061405429*t**5)*Exp(-x**2))*(x>0)-\
		(1-(0.254829592*t-0.284496736*t**2+1.421413741*t**3-1.453152027*t**4+1.061405429*t**5)*Exp(-x**2))*(x<0)
        Pfail=0.5+0.5*Erf_1 #static probability
        arcpy.AddMessage("Probability of FoS < 1 Created")

        arcpy.AddMessage("Calculating Reliability Index")
        
        RI=(FS-1)/FS_SD #static reliability
		
        Pfail.save("Pfail")
        RI.save("RI")
        arcpy.AddMessage("Reliability Index Created")

        arcpy.AddMessage("Calculating Critical Newmark acceleration and SD")
		
        aNcrit=(FS-1)*Sin(beta) # seismic_a mean
        aNcrit_SD=SquareRoot((FS-1)**2*(Cos(beta))**2*beta_SD**2+(Sin(beta))**2*FS_SD**2)
        
        aNcrit.save("aNcrit") # seismic_a sd
        aNcrit_SD.save("aNcrit_SD")
        arcpy.AddMessage("Critical Newmark acceleration and SD created")

        arcpy.AddMessage("Calculating Probability of exceeding critical Newmark acceleration")
		
		# Following Abramowitz and Stegun (1964), equation 7.1.26 (pp 299), with an upper error bound of 1.5x10^-7:
        x=(SeisA-aNcrit)/(aNcrit_SD*SquareRoot(2))
        t=1/(1+0.3275911*Abs(x))
        Erf_1=(1-(0.254829592*t-0.284496736*t**2+1.421413741*t**3-1.453152027*t**4+1.061405429*t**5)*Exp(-x**2))*(x>0)-\
		(1-(0.254829592*t-0.284496736*t**2+1.421413741*t**3-1.453152027*t**4+1.061405429*t**5)*Exp(-x**2))*(x<0)
        PSeisA=0.5+0.5*Erf_1 #seismic_a probability
        PSeisA.save("PSeisA")
        arcpy.AddMessage("Probability of exceeding critical Newmark acceleration Created")
		
        arcpy.AddMessage("Calculating Reliability for critical Newmark acceleration")
		
        RISeis=(aNcrit-SeisA)/aNcrit_SD #seismic_a reliability
        RISeis.save("RISeis")
        arcpy.AddMessage("Reliability for critical Newmark acceleration Created")
		
		
        arcpy.AddMessage("Deleting temporary files")
		
        if SlopeMeth=="second-order accurate central difference approximation (PISA-m method)":

            arcpy.Delete_management("DEM_N")
            arcpy.Delete_management("DEM_S")
            arcpy.Delete_management("DEM_E")
            arcpy.Delete_management("DEM_W")
			
        elif SlopeMeth=="Eight neighbour pixel method (ArcMap method)":
			
            arcpy.Delete_management("DEM_z1")
            arcpy.Delete_management("DEM_z2")
            arcpy.Delete_management("DEM_z3")
            arcpy.Delete_management("DEM_z4")
            arcpy.Delete_management("DEM_z5")
            arcpy.Delete_management("DEM_z6")
            arcpy.Delete_management("DEM_z7")
            arcpy.Delete_management("DEM_z8")
					
        arcpy.AddMessage("*****  All calculations successfuly finished  *****")
		
		
        arcpy.AddMessage("*****  All calculations successfuly finished  *****")

		
        return
		

		
		
		
class GIS_TISSA_GUI(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "GIS_TISSA_GUI"
        self.description = "GIS_TISSA reading inputs through the GUI"
        # self.description = "Some description..."
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""

		# DEM input
        param0 = arcpy.Parameter(
		    name="Digital Elevation Model",
            displayName="Digital Elevation Model",
            datatype="DERasterDataset",
            parameterType="Required",
            direction="Input"
        )
        param0.value="dem.tif"
		
		# DEM standard deviation input
        param1 = arcpy.Parameter(
		    name="DEM elevation standard deviation",
            displayName="DEM elevation sTandard deviation",
            datatype="GPDouble",
            parameterType="Required",
            direction="Input"
        )
        param1.value=2

		
		# Phi inputs
        param2 = arcpy.Parameter(
		    name="phi_dist",
            displayName="Distribution type for the internal friction angle",
            # datatype="DEFeatureClass",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param2.filter.type="ValueList"
        param2.filter.list=["Normal", "Uniform", "Gumbel", "Triangular","Beta-PERT"]
        param2.value="Normal"

        param3 = arcpy.Parameter(
		    name="Phi_dist_param_1",
            displayName="Parameter 1 of the distribution of the Internal friction angle",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Required",
            direction="Input"
        )
        param3.value="phi_p1.tif"

        param4 = arcpy.Parameter(
		    name="Phi_dist_param_2",
            displayName="Parameter 2 of the distribution of the Internal friction angle",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )
        param4.value="phi_p2.tif"

        param5 = arcpy.Parameter(
		    name="Phi_dist_param_3",
            displayName="Parameter 3 of the distribution of the Internal friction angle",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )

		
		# C inputs
        param6 = arcpy.Parameter(
		    name="c_dist",
            displayName="Distribution type for the soil cohesion",
            # datatype="DEFeatureClass",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param6.filter.type="ValueList"
        param6.filter.list=["Normal", "Uniform", "Gumbel", "Triangular","Beta-PERT"]
        param6.value="Normal"

        param7 = arcpy.Parameter(
		    name="c_dist_param_1",
            displayName="Parameter 1 of the distribution of the soil cohesion",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Required",
            direction="Input"
        )
        param7.value="c_p1.tif"

        param8 = arcpy.Parameter(
		    name="c_dist_param_2",
            displayName="Parameter 2 of the distribution of the soil cohesion",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )
        param8.value="c_p2.tif"

        param9 = arcpy.Parameter(
		    name="c_dist_param_3",
            displayName="Parameter 3 of the distribution of the soil cohesion",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )		
		
		

		# gm inputs
        param10 = arcpy.Parameter(
		    name="gm_dist",
            displayName="Distribution type for the moist soil weight",
            # datatype="DEFeatureClass",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param10.filter.type="ValueList"
        param10.filter.list=["Normal", "Uniform", "Gumbel", "Triangular","Beta-PERT"]
        param10.value="Normal"

        param11 = arcpy.Parameter(
		    name="gm_dist_param_1",
            displayName="Parameter 1 of the distribution of the moist soil weight",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Required",
            direction="Input"
        )
        param11.value="muw_p1.tif"

        param12 = arcpy.Parameter(
		    name="gm_dist_param_2",
            displayName="Parameter 2 of the distribution of the moist soil weight",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )
        param12.value="muw_p2.tif"

        param13 = arcpy.Parameter(
		    name="gm_dist_param_3",
            displayName="Parameter 3 of the distribution of the moist soil weight",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )		

		
		
		
		
		# gs inputs
        param14 = arcpy.Parameter(
		    name="gs_dist",
            displayName="Distribution type for the saturated soil weight",
            # datatype="DEFeatureClass",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param14.filter.type="ValueList"
        param14.filter.list=["Normal", "Uniform", "Gumbel", "Triangular","Beta-PERT"]
        param14.value="Normal"

        param15 = arcpy.Parameter(
		    name="gs_dist_param_1",
            displayName="Parameter 1 of the distribution of the saturated soil weight",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Required",
            direction="Input"
        )
        param15.value="suw_p1.tif"

        param16 = arcpy.Parameter(
		    name="gs_dist_param_2",
            displayName="Parameter 2 of the distribution of the saturated soil weight",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )
        param16.value="suw_p2.tif"

        param17 = arcpy.Parameter(
		    name="gs_dist_param_3",
            displayName="Parameter 3 of the distribution of the saturated soil weight",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )		




		# D inputs
        param18 = arcpy.Parameter(
		    name="D_dist",
            displayName="Distribution type for the soil depth",
            # datatype="DEFeatureClass",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param18.filter.type="ValueList"
        param18.filter.list=["Normal", "Uniform", "Gumbel", "Triangular","Beta-PERT"]
        param18.value="Normal"

        param19 = arcpy.Parameter(
		    name="D_dist_param_1",
            displayName="Parameter 1 of the distribution of the soil depth",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Required",
            direction="Input"
        )
        param19.value="D_p1.tif"

        param20 = arcpy.Parameter(
		    name="D_dist_param_2",
            displayName="Parameter 2 of the distribution of the soil depth",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )
        param20.value="D_p2.tif"

        param21 = arcpy.Parameter(
		    name="D_dist_param_3",
            displayName="Parameter 3 of the distribution of the soil depth",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )		

		
		
		
		
		# Hw inputs
        param22 = arcpy.Parameter(
		    name="Hw_dist",
            displayName="Distribution type for the soil phreatic ratio",
            # datatype="DEFeatureClass",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param22.filter.type="ValueList"
        param22.filter.list=["Normal", "Uniform", "Gumbel", "Triangular","Beta-PERT"]
        param22.value="Normal"

        param23 = arcpy.Parameter(
		    name="Hw_dist_param_1",
            displayName="Parameter 1 of the distribution of the soil phreatic ratio",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Required",
            direction="Input"
        )
        param23.value="Hw_p1.tif"

        param24 = arcpy.Parameter(
		    name="Hw_dist_param_2",
            displayName="Parameter 2 of the distribution of the soil phreatic ratio",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )
        param24.value="Hw_p2.tif"

        param25 = arcpy.Parameter(
		    name="Hw_dist_param_3",
            displayName="Parameter 3 of the distribution of the soil phreatic ratio",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )		



		# Cr inputs
        param26 = arcpy.Parameter(
		    name="cr_dist",
            displayName="Distribution type for the root cohesion",
            # datatype="DEFeatureClass",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param26.filter.type="ValueList"
        param26.filter.list=["Normal", "Uniform", "Gumbel", "Triangular","Beta-PERT"]
        param26.value="Normal"

        param27 = arcpy.Parameter(
		    name="cr_dist_param_1",
            displayName="Parameter 1 of the distribution of the root cohesion",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Required",
            direction="Input"
        )
        param27.value="Cr_p1.tif"

        param28 = arcpy.Parameter(
		    name="cr_dist_param_2",
            displayName="Parameter 2 of the distribution of the root cohesion",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )
        param28.value="Cr_p2.tif"

        param29 = arcpy.Parameter(
		    name="cr_dist_param_3",
            displayName="Parameter 3 of the distribution of the root cohesion",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )		





		# qt inputs
        param30 = arcpy.Parameter(
		    name="qt_dist",
            displayName="Distribution type for the vegetation surcharge",
            # datatype="DEFeatureClass",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param30.filter.type="ValueList"
        param30.filter.list=["Normal", "Uniform", "Gumbel", "Triangular","Beta-PERT"]
        param30.value="Normal"

        param31 = arcpy.Parameter(
		    name="qt_dist_param_1",
            displayName="Parameter 1 of the distribution of the vegetation surcharge",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Required",
            direction="Input"
        )
        param31.value="qt_p1.tif"

        param32 = arcpy.Parameter(
		    name="qt_dist_param_2",
            displayName="Parameter 2 of the distribution of the vegetation surcharge",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )
        param32.value="qt_p2.tif"

        param33 = arcpy.Parameter(
		    name="qt_dist_param_3",
            displayName="Parameter 3 of the distribution of the vegetation surcharge",
            datatype=["DERasterDataset","GPDouble"],
            parameterType="Optional",
            direction="Input"
        )		


		
		
        param34 = arcpy.Parameter(
		    name="Minimum slope value",
            displayName="Minimum Slope Value",
            datatype="GPDouble",
            parameterType="Required",
            direction="Input"
        )
        param34.value=5

        param35 = arcpy.Parameter(
		    name="Output geodatabase",
            displayName="Output Geodataset",
            datatype="DEWorkspace",
            parameterType="Required",
            direction="Input"
        )

        param36 = arcpy.Parameter(
		    name="Expected seismic ground acceleration",
            displayName="Expected seismic ground acceleration",
            datatype="GPDouble",
            parameterType="Required",
            direction="Input"
        )
        param36.value=0.39
		
		
        param39 = arcpy.Parameter(
		    name="Slope_method",
            displayName="Slope calculation method",
            datatype="String",
            parameterType="Optional",
            direction="Input"
        )
        param39.filter.type="ValueList"
        param39.filter.list=["second-order accurate central difference approximation (PISA-m method)", "Eight neighbour pixel method (ArcMap method)"]
        param39.value="Eight neighbour pixel method (ArcMap method)"
		
		
        # params = [param0, param1, param2, param3, param4, param5, param6, param7, param8, param9,
                  # param10]
        params = [param0, param1, param2, param3, param4, param5, param6, param7, param8, param9,
		param10, param11, param12, param13, param14, param15, param16, param17, param18, param19,
		param20, param21, param22, param23, param24, param25, param26, param27, param28, param29,
		param30, param31, param32, param33, param34, param35, param36, param39]
				  
        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def execute(self, parameters, messages):
        """The source code of the tool."""
        arcpy.CheckOutExtension("Spatial")
        arcpy.CheckOutExtension("3D")

        arcpy.AddMessage("Assigning inputs")
		
		# Assigning inputs
        DEM_1 = parameters[0].valueAsText
        DEM_SD = parameters[1].valueAsText
        phi_dist = parameters[2].valueAsText
        phi_p1 = parameters[3].valueAsText
        phi_p2 = parameters[4].valueAsText
        phi_p3 = parameters[5].valueAsText
        c_dist = parameters[6].valueAsText
        c_p1 = parameters[7].valueAsText
        c_p2 = parameters[8].valueAsText
        c_p3 = parameters[9].valueAsText
        gm_dist = parameters[10].valueAsText
        gm_p1 = parameters[11].valueAsText
        gm_p2 = parameters[12].valueAsText
        gm_p3 = parameters[13].valueAsText
        gs_dist = parameters[14].valueAsText
        gs_p1 = parameters[15].valueAsText
        gs_p2 = parameters[16].valueAsText
        gs_p3 = parameters[17].valueAsText
        D_dist = parameters[18].valueAsText
        D_p1 = parameters[19].valueAsText
        D_p2 = parameters[20].valueAsText
        D_p3 = parameters[21].valueAsText
        Hw_dist = parameters[22].valueAsText
        Hw_p1 = parameters[23].valueAsText
        Hw_p2 = parameters[24].valueAsText
        Hw_p3 = parameters[25].valueAsText
        Cr_dist = parameters[26].valueAsText
        Cr_p1 = parameters[27].valueAsText
        Cr_p2 = parameters[28].valueAsText
        Cr_p3 = parameters[29].valueAsText
        qt_dist = parameters[30].valueAsText
        qt_p1 = parameters[31].valueAsText
        qt_p2 = parameters[32].valueAsText
        qt_p3 = parameters[33].valueAsText
        Minimum_slope_value = parameters[34].valueAsText
        wd=parameters[35].valueAsText
        SeisA0 = parameters[36].valueAsText
        SlopeMeth = parameters[37].valueAsText

        DEM_2=Raster(DEM_1)
		
		# Define workspace and overwritting permission
        arcpy.env.workspace=wd
        arcpy.env.scratchWorkspace=wd
        arcpy.env.overwriteOutput = True

		# Extract the pixel size and raste extent from the DEM
        cellsz=arcpy.GetRasterProperties_management(DEM_2,"CELLSIZEX")
        x1=arcpy.GetRasterProperties_management(DEM_1,"LEFT")
        x2=arcpy.GetRasterProperties_management(DEM_1,"RIGHT")
        y1=arcpy.GetRasterProperties_management(DEM_1,"BOTTOM")
        y2=arcpy.GetRasterProperties_management(DEM_1,"TOP")

		#Create constant rasters 
        gw=CreateConstantRaster(9810,"Float",cellsz,Extent(x1, y1, x2, y2))
        SeisA=CreateConstantRaster(SeisA0,"Float",cellsz,Extent(x1, y1, x2, y2))
		
		
        if SlopeMeth=="second-order accurate central difference approximation (PISA-m method)":

            arcpy.AddMessage("Slope calculation compatible with PISA-m")
		
		    # Slope calculation compatible with PISA-m
            cellsz2=cellsz.getOutput(0)
            cellsz4=float(cellsz2)
            nshift="-"+cellsz.getOutput(0)
		
            arcpy.Shift_management(DEM_2,"DEM_N",0,cellsz)
            arcpy.Shift_management(DEM_2,"DEM_S",0,nshift)
            arcpy.Shift_management(DEM_2,"DEM_E",cellsz,0)
            arcpy.Shift_management(DEM_2,"DEM_W",nshift,0)
		
            DEM_N2=Raster("DEM_N")
            DEM_S2=Raster("DEM_S")
            DEM_E2=Raster("DEM_E")
            DEM_W2=Raster("DEM_W")

            beta=ATan(SquareRoot((DEM_N2-DEM_S2)**2+(DEM_E2-DEM_W2)**2)/(2*cellsz4))
            beta_SD=SquareRoot((8*cellsz4**2*float(DEM_SD)**2)/((4*cellsz4**2+(DEM_N2-DEM_S2)**2+(DEM_E2-DEM_W2)**2)**2))
			
			
            # arcpy.Delete_management("DEM_N")
            # arcpy.Delete_management("DEM_S")
            # arcpy.Delete_management("DEM_E")
            # arcpy.Delete_management("DEM_W")

			
        elif SlopeMeth=="Eight neighbour pixel method (ArcMap method)":
		
            arcpy.AddMessage("Standard ArcMap slope calculation")
			
            beta=Slope(DEM_2,"DEGREE")*3.1415926/180		
			
		    # Slope calculation compatible with PISA-m
            cellsz2=cellsz.getOutput(0)
            cellsz4=float(cellsz2)
            nshift="-"+cellsz.getOutput(0)
		
            arcpy.Shift_management(DEM_2,"DEM_z1",cellsz,cellsz)
            arcpy.Shift_management(DEM_2,"DEM_z2",cellsz,0)
            arcpy.Shift_management(DEM_2,"DEM_z3",cellsz,nshift)
            arcpy.Shift_management(DEM_2,"DEM_z4",nshift,cellsz)
            arcpy.Shift_management(DEM_2,"DEM_z5",nshift,0)
            arcpy.Shift_management(DEM_2,"DEM_z6",nshift,nshift)
            arcpy.Shift_management(DEM_2,"DEM_z7",0,nshift)
            arcpy.Shift_management(DEM_2,"DEM_z8",0,cellsz)
		
            DEM_Z1=Raster("DEM_z1")
            DEM_Z2=Raster("DEM_z2")
            DEM_Z3=Raster("DEM_z3")
            DEM_Z4=Raster("DEM_z4")
            DEM_Z5=Raster("DEM_z5")
            DEM_Z6=Raster("DEM_z6")
            DEM_Z7=Raster("DEM_z7")
            DEM_Z8=Raster("DEM_z8")
			
            ADEM=((DEM_Z1+2*DEM_Z2+DEM_Z3)-(DEM_Z4+2*DEM_Z5+DEM_Z6))/(8*cellsz4)
            BDEM=((DEM_Z6+2*DEM_Z7+DEM_Z3)-(DEM_Z4+2*DEM_Z8+DEM_Z1))/(8*cellsz4)
			
            beta_SD=SquareRoot(3*float(DEM_SD)**2/(8*cellsz4**2)*(ADEM**2+BDEM**2)/((1+ADEM**2+BDEM**2)*SquareRoot(ADEM**2+BDEM**2)))		

			
            # arcpy.Delete_management("DEM_z1")
            # arcpy.Delete_management("DEM_z2")
            # arcpy.Delete_management("DEM_z3")
            # arcpy.Delete_management("DEM_z4")
            # arcpy.Delete_management("DEM_z5")
            # arcpy.Delete_management("DEM_z6")
            # arcpy.Delete_management("DEM_z7")
            # arcpy.Delete_management("DEM_z8")

			

		# Convert phi inputs to rasters
        arcpy.AddMessage("Convert phi inputs to rasters")
        if phi_dist=="Triangular" or phi_dist=="Beta-PERT":
            if not arcpy.Exists(parameters[5].value):
                phi_p3=CreateConstantRaster(phi_p3,"Float",cellsz,Extent(x1, y1, x2, y2))
        if not arcpy.Exists(parameters[3].value):
            phi_p1=CreateConstantRaster(phi_p1,"Float",cellsz,Extent(x1, y1, x2, y2))		
        if not arcpy.Exists(parameters[4].value):
            phi_p2=CreateConstantRaster(phi_p2,"Float",cellsz,Extent(x1, y1, x2, y2))
			
		# Calculate mean and standard deviation for phi
        arcpy.AddMessage("Calculate mean and standard deviation for phi")
        if phi_dist=="Normal":
            phi4=phi_p1
            phi4_SD=phi_p2
        elif phi_dist=="Uniform":
            phi4=(phi_p1+phi_p2)/2
            phi4_SD=(phi_p1-phi_p2)**2/2
        elif phi_dist=="Gumbel":
            phi4=phi_p1+0.577216*phi_p2
            phi4_SD=(3.1415926**2)*phi_p2**2/6
        elif phi_dist=="Triangular":
            phi4=(phi_p1+phi_p2+phi_p3)/3
            phi4_SD=(phi_p1*(phi_p1-phi_p2)+phi_p2*(phi_p2-phi_p3)+phi_p3*(phi_p3-phi_p1))/18
        elif phi_dist=="Beta-PERT":
            phi4=(phi_p1+4*phi_p2+phi_p3)/6
            phi4_SD=(phi_p1-phi_p3)**2/36
			
		# Convert Cs inputs to rasters
        arcpy.AddMessage("Convert Cs inputs to rasters")
        if c_dist=="Triangular" or c_dist=="Beta-PERT":
            if not arcpy.Exists(parameters[9].value):
                c_p3=CreateConstantRaster(c_p3,"Float",cellsz,Extent(x1, y1, x2, y2))
        if not arcpy.Exists(parameters[7].value):
            c_p1=CreateConstantRaster(c_p1,"Float",cellsz,Extent(x1, y1, x2, y2))		
        if not arcpy.Exists(parameters[8].value):
            c_p2=CreateConstantRaster(c_p2,"Float",cellsz,Extent(x1, y1, x2, y2))
			
		# Calculate mean and standard deviation for Cs
        arcpy.AddMessage("Calculate mean and standard deviation for Cs")
        if c_dist=="Normal":
            Cs=c_p1
            Cs_SD=c_p2
        elif c_dist=="Uniform":
            Cs=(c_p1+c_p2)/2
            Cs_SD=(c_p1-c_p2)**2/2
        elif c_dist=="Gumbel":
            Cs=c_p1+0.577216*c_p2
            Cs_SD=(3.1415926**2)*c_p2**2/6
        elif c_dist=="Triangular":
            Cs=(c_p1+c_p2+c_p3)/3
            Cs_SD=(c_p1*(c_p1-c_p2)+c_p2*(c_p2-c_p3)+c_p3*(c_p3-c_p1))/18
        elif c_dist=="Beta-PERT":
            Cs=(c_p1+4*c_p2+c_p3)/6
            Cs_SD=(c_p1-c_p3)**2/36

		# Convert gm inputs to rasters
        arcpy.AddMessage("Convert gm inputs to rasters")
        if gm_dist=="Triangular" or gm_dist=="Beta-PERT":
            if not arcpy.Exists(parameters[13].value):
                gm_p3=CreateConstantRaster(gm_p3,"Float",cellsz,Extent(x1, y1, x2, y2))
        if not arcpy.Exists(parameters[11].value):
            gm_p1=CreateConstantRaster(gm_p1,"Float",cellsz,Extent(x1, y1, x2, y2))		
        if not arcpy.Exists(parameters[12].value):
            gm_p2=CreateConstantRaster(gm_p2,"Float",cellsz,Extent(x1, y1, x2, y2))
			
		# Calculate mean and standard deviation for gm
        arcpy.AddMessage("Calculate mean and standard deviation for gm")
        if gm_dist=="Normal":
            gm=gm_p1
            gm_SD=gm_p2
        elif gm_dist=="Uniform":
            gm=(gm_p1+gm_p2)/2
            gm_SD=(gm_p1-gm_p2)**2/2
        elif gm_dist=="Gumbel":
            gm=gm_p1+0.577216*gm_p2
            gm_SD=(3.1415926**2)*gm_p2**2/6
        elif gm_dist=="Triangular":
            gm=(gm_p1+gm_p2+gm_p3)/3
            gm_SD=(gm_p1*(gm_p1-gm_p2)+gm_p2*(gm_p2-gm_p3)+gm_p3*(gm_p3-gm_p1))/18
        elif gm_dist=="Beta-PERT":
            gm=(gm_p1+4*gm_p2+gm_p3)/6
            gm_SD=(gm_p1-gm_p3)**2/36


		# Convert gs inputs to rasters
        arcpy.AddMessage("Convert gs inputs to rasters")
        if gs_dist=="Triangular" or gs_dist=="Beta-PERT":
            if not arcpy.Exists(parameters[17].value):
                gs_p3=CreateConstantRaster(gs_p3,"Float",cellsz,Extent(x1, y1, x2, y2))
        if not arcpy.Exists(parameters[15].value):
            gs_p1=CreateConstantRaster(gs_p1,"Float",cellsz,Extent(x1, y1, x2, y2))		
        if not arcpy.Exists(parameters[16].value):
            gs_p2=CreateConstantRaster(gs_p2,"Float",cellsz,Extent(x1, y1, x2, y2))
			
		# Calculate mean and standard deviation for gs
        arcpy.AddMessage("Calculate mean and standard deviation for gs")
        if gs_dist=="Normal":
            gs=gs_p1
            gs_SD=gs_p2
        elif gs_dist=="Uniform":
            gs=(gs_p1+gs_p2)/2
            gs_SD=(gs_p1-gs_p2)**2/2
        elif gs_dist=="Gumbel":
            gs=gs_p1+0.577216*gs_p2
            gs_SD=(3.1415926**2)*gs_p2**2/6
        elif gs_dist=="Triangular":
            gs=(gs_p1+gs_p2+gs_p3)/3
            gs_SD=(gs_p1*(gs_p1-gs_p2)+gs_p2*(gs_p2-gs_p3)+gs_p3*(gs_p3-gs_p1))/18
        elif gs_dist=="Beta-PERT":
            gs=(gs_p1+4*gs_p2+gs_p3)/6
            gs_SD=(gs_p1-gs_p3)**2/36


		# Convert D inputs to rasters
        arcpy.AddMessage("Convert D inputs to rasters")
        if D_dist=="Triangular" or D_dist=="Beta-PERT":
            if not arcpy.Exists(parameters[21].value):
                D_p3=CreateConstantRaster(D_p3,"Float",cellsz,Extent(x1, y1, x2, y2))
        if not arcpy.Exists(parameters[19].value):
            D_p1=CreateConstantRaster(D_p1,"Float",cellsz,Extent(x1, y1, x2, y2))		
        if not arcpy.Exists(parameters[20].value):
            D_p2=CreateConstantRaster(D_p2,"Float",cellsz,Extent(x1, y1, x2, y2))
			
		# Calculate mean and standard deviation for D
        arcpy.AddMessage("Calculate mean and standard deviation for D")
        if D_dist=="Normal":
            D=D_p1
            D_SD=D_p2
        elif D_dist=="Uniform":
            D=(D_p1+D_p2)/2
            D_SD=(D_p1-D_p2)**2/2
        elif D_dist=="Gumbel":
            D=D_p1+0.577216*D_p2
            D_SD=(3.1415926**2)*D_p2**2/6
        elif D_dist=="Triangular":
            D=(D_p1+D_p2+D_p3)/3
            D_SD=(D_p1*(D_p1-D_p2)+D_p2*(D_p2-D_p3)+D_p3*(D_p3-D_p1))/18
        elif D_dist=="Beta-PERT":
            D=(D_p1+4*D_p2+D_p3)/6
            D_SD=(D_p1-D_p3)**2/36


		# Convert Hw inputs to rasters
        arcpy.AddMessage("Convert Hw inputs to rasters")
        if Hw_dist=="Triangular" or Hw_dist=="Beta-PERT":
            if not arcpy.Exists(parameters[25].value):
                Hw_p3=CreateConstantRaster(Hw_p3,"Float",cellsz,Extent(x1, y1, x2, y2))
        if not arcpy.Exists(parameters[23].value):
            Hw_p1=CreateConstantRaster(Hw_p1,"Float",cellsz,Extent(x1, y1, x2, y2))		
        if not arcpy.Exists(parameters[24].value):
            Hw_p2=CreateConstantRaster(Hw_p2,"Float",cellsz,Extent(x1, y1, x2, y2))
			
		# Calculate mean and standard deviation for Hw
        arcpy.AddMessage("Calculate mean and standard deviation for Hw")
        if Hw_dist=="Normal":
            Hw=Hw_p1
            Hw_SD=Hw_p2
        elif Hw_dist=="Uniform":
            Hw=(Hw_p1+Hw_p2)/2
            Hw_SD=(Hw_p1-Hw_p2)**2/2
        elif Hw_dist=="Gumbel":
            Hw=Hw_p1+0.577216*Hw_p2
            Hw_SD=(3.1415926**2)*Hw_p2**2/6
        elif Hw_dist=="Triangular":
            Hw=(Hw_p1+Hw_p2+Hw_p3)/3
            Hw_SD=(Hw_p1*(Hw_p1-Hw_p2)+Hw_p2*(Hw_p2-Hw_p3)+Hw_p3*(Hw_p3-Hw_p1))/18
        elif Hw_dist=="Beta-PERT":
            Hw=(Hw_p1+4*Hw_p2+Hw_p3)/6
            Hw_SD=(Hw_p1-Hw_p3)**2/36


		# Convert Cr inputs to rasters
        arcpy.AddMessage("Convert Cr inputs to rasters")
        if Cr_dist=="Triangular" or Cr_dist=="Beta-PERT":
            if not arcpy.Exists(parameters[29].value):
                Cr_p3=CreateConstantRaster(Cr_p3,"Float",cellsz,Extent(x1, y1, x2, y2))
        if not arcpy.Exists(parameters[27].value):
            Cr_p1=CreateConstantRaster(Cr_p1,"Float",cellsz,Extent(x1, y1, x2, y2))		
        if not arcpy.Exists(parameters[28].value):
            Cr_p2=CreateConstantRaster(Cr_p2,"Float",cellsz,Extent(x1, y1, x2, y2))
			
		# Calculate mean and standard deviation for Cr
        arcpy.AddMessage("Calculate mean and standard deviation for Cr")
        if Cr_dist=="Normal":
            Cr=Cr_p1
            Cr_SD=Cr_p2
        elif Cr_dist=="Uniform":
            Cr=(Cr_p1+Cr_p2)/2
            Cr_SD=(Cr_p1-Cr_p2)**2/2
        elif Cr_dist=="Gumbel":
            Cr=Cr_p1+0.577216*Cr_p2
            Cr_SD=(3.1415926**2)*Cr_p2**2/6
        elif Cr_dist=="Triangular":
            Cr=(Cr_p1+Cr_p2+Cr_p3)/3
            Cr_SD=(Cr_p1*(Cr_p1-Cr_p2)+Cr_p2*(Cr_p2-Cr_p3)+Cr_p3*(Cr_p3-Cr_p1))/18
        elif Cr_dist=="Beta-PERT":
            Cr=(Cr_p1+4*Cr_p2+Cr_p3)/6
            Cr_SD=(Cr_p1-Cr_p3)**2/36


		# Convert qt inputs to rasters
        arcpy.AddMessage("Convert qt inputs to rasters")
        if qt_dist=="Triangular" or qt_dist=="Beta-PERT":
            if not arcpy.Exists(parameters[33].value):
                qt_p3=CreateConstantRaster(qt_p3,"Float",cellsz,Extent(x1, y1, x2, y2))
        if not arcpy.Exists(parameters[31].value):
            qt_p1=CreateConstantRaster(qt_p1,"Float",cellsz,Extent(x1, y1, x2, y2))		
        if not arcpy.Exists(parameters[32].value):
            qt_p2=CreateConstantRaster(qt_p2,"Float",cellsz,Extent(x1, y1, x2, y2))
			
		# Calculate mean and standard deviation for qt
        arcpy.AddMessage("Calculate mean and standard deviation for qt")
        if qt_dist=="Normal":
            qt=qt_p1
            qt_SD=qt_p2
        elif qt_dist=="Uniform":
            qt=(qt_p1+qt_p2)/2
            qt_SD=(qt_p1-qt_p2)**2/2
        elif qt_dist=="Gumbel":
            qt=qt_p1+0.577216*qt_p2
            qt_SD=(3.1415926**2)*qt_p2**2/6
        elif qt_dist=="Triangular":
            qt=(qt_p1+qt_p2+qt_p3)/3
            qt_SD=(qt_p1*(qt_p1-qt_p2)+qt_p2*(qt_p2-qt_p3)+qt_p3*(qt_p3-qt_p1))/18
        elif qt_dist=="Beta-PERT":
            qt=(qt_p1+4*qt_p2+qt_p3)/6
            qt_SD=(qt_p1-qt_p3)**2/36


				
        phi=phi4*CreateConstantRaster(3.1415926/180,"Float",cellsz,Extent(x1, y1, x2, y2))
        # phi.save("phi")
        phi_SD=phi4_SD*CreateConstantRaster(3.1415926/180,"Float",cellsz,Extent(x1, y1, x2, y2))
        # phi_SD.save("phi_SD")

        arcpy.AddMessage("Calculating Mean Factor of Safety ")
		
        Cr=Raster(Cr)
        Cs=Raster(Cs)
        qt=Raster(qt)
        gm=Raster(gm)
        D=Raster(D)
        gs=Raster(gs)
        # gw=Raster(gw)
        Hw=Raster(Hw)
        # beta=Raster(beta)
        # phi=Raster(phi)
		
        FS=(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))
        FS.save("FoS_Mean") #static mean
        arcpy.AddMessage("Mean Factor of Safety Calculated")

        arcpy.AddMessage("Calculating Standard Deviation of Factor of Safety")
		
		
		
        Dcr=1/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))
        Dcs=1/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))
        Dqt=((Cos(beta))**2*Tan(phi)*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-Sin(beta)*Cos(beta)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        Dgm=(((Cos(beta))**2*Tan(phi)*(D-Hw*D))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-Sin(beta)*Cos(beta)*(D-Hw*D)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        Dgs=(((Cos(beta))**2*Tan(phi)*(Hw*D))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-Sin(beta)*Cos(beta)*(Hw*D)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        DD=(((Cos(beta))**2*Tan(phi)*(gm+(gs-gw-gm)*Hw))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-Sin(beta)*Cos(beta)*(gm+(gs-gm)*Hw)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        DHw=(((gs-gw-gm)*D*(Cos(beta))**2*Tan(phi))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-(gs-gm)*D*Sin(beta)*Cos(beta)*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        Dbeta=((2*Cos(beta)*(-Sin(beta))*(qt+gm*D+(gs-gw-gm)*Hw*D)*Tan(phi))*((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))-(qt+gm*D+(gs-gm)*Hw*D)*((Cos(beta))**2+(-Sin(beta))*Sin(beta))*(Cr+Cs+(qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*Tan(phi)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))**2
        Dphi=((qt+gm*D+(gs-gw-gm)*Hw*D)*(Cos(beta))**2*(1/(Cos(phi)**2)))/((qt+gm*D+(gs-gm)*Hw*D)*Sin(beta)*Cos(beta))
		
        Cr_SD=Raster(Cr_SD)
        Cs_SD=Raster(Cs_SD)
        qt_SD=Raster(qt_SD)
        gm_SD=Raster(gm_SD)
        D_SD=Raster(D_SD)
        gs_SD=Raster(gs_SD)
        Hw_SD=Raster(Hw_SD)
        # beta_SD=Raster(beta_SD)
        # phi_SD=Raster(phi_SD)
		
		
        FS_SD=SquareRoot(Dcr**2*Cr_SD**2+Dcs**2*Cs_SD**2+Dqt**2*qt_SD**2+Dgm**2*gm_SD**2+Dgs**2*gs_SD**2+DD**2*D_SD**2+DHw**2*Hw_SD**2+Dbeta**2*beta_SD**2+Dphi**2*phi_SD**2)
        FS_SD.save("FoS_SD")	# static sd	
        arcpy.AddMessage("Standard Deviation of Factor of Safety Created")

        arcpy.AddMessage("Calculating Probability of FoS < 1")
		
		# Following Abramowitz and Stegun (1964), equation 7.1.26 (pp 299), with an upper error bound of 1.5x10^-7:
        x=-1*(Ln(FS)-0.5*Ln((FS_SD/FS)**2+1))/(SquareRoot(2)*SquareRoot(Ln((FS_SD/FS)**2+1)))
        t=1/(1+0.3275911*Abs(x))
        Erf_1=(1-(0.254829592*t-0.284496736*t**2+1.421413741*t**3-1.453152027*t**4+1.061405429*t**5)*Exp(-x**2))*(x>0)-\
		(1-(0.254829592*t-0.284496736*t**2+1.421413741*t**3-1.453152027*t**4+1.061405429*t**5)*Exp(-x**2))*(x<0)
        Pfail=0.5+0.5*Erf_1 #static probability
        arcpy.AddMessage("Probability of FoS < 1 Created")

        arcpy.AddMessage("Calculating Reliability Index")
        
        RI=(FS-1)/FS_SD #static reliability
		
        Pfail.save("Pfail")
        RI.save("RI")
        arcpy.AddMessage("Reliability Index Created")

        arcpy.AddMessage("Calculating Critical Newmark acceleration and SD")
		
        aNcrit=(FS-1)*Sin(beta) # seismic_a mean
        aNcrit_SD=SquareRoot((FS-1)**2*(Cos(beta))**2*beta_SD**2+(Sin(beta))**2*FS_SD**2)
        
        aNcrit.save("aNcrit") # seismic_a sd
        aNcrit_SD.save("aNcrit_SD")
        arcpy.AddMessage("Critical Newmark acceleration and SD created")

        arcpy.AddMessage("Calculating Probability of exceeding critical Newmark acceleration")
		
		# Following Abramowitz and Stegun (1964), equation 7.1.26 (pp 299), with an upper error bound of 1.5x10^-7:
        x=(SeisA-aNcrit)/(aNcrit_SD*SquareRoot(2))
        t=1/(1+0.3275911*Abs(x))
        Erf_1=(1-(0.254829592*t-0.284496736*t**2+1.421413741*t**3-1.453152027*t**4+1.061405429*t**5)*Exp(-x**2))*(x>0)-\
		(1-(0.254829592*t-0.284496736*t**2+1.421413741*t**3-1.453152027*t**4+1.061405429*t**5)*Exp(-x**2))*(x<0)
        PSeisA=0.5+0.5*Erf_1 #seismic_a probability
        PSeisA.save("PSeisA")
        arcpy.AddMessage("Probability of exceeding critical Newmark acceleration Created")
		
        arcpy.AddMessage("Calculating Reliability for critical Newmark acceleration")
		
        RISeis=(aNcrit-SeisA)/aNcrit_SD #seismic_a reliability
        RISeis.save("RISeis")
        arcpy.AddMessage("Reliability for critical Newmark acceleration Created")
		
		
		
        arcpy.AddMessage("Deleting temporary files")
		
        if SlopeMeth=="second-order accurate central difference approximation (PISA-m method)":

            arcpy.Delete_management("DEM_N")
            arcpy.Delete_management("DEM_S")
            arcpy.Delete_management("DEM_E")
            arcpy.Delete_management("DEM_W")
			
        elif SlopeMeth=="Eight neighbour pixel method (ArcMap method)":
			
            arcpy.Delete_management("DEM_z1")
            arcpy.Delete_management("DEM_z2")
            arcpy.Delete_management("DEM_z3")
            arcpy.Delete_management("DEM_z4")
            arcpy.Delete_management("DEM_z5")
            arcpy.Delete_management("DEM_z6")
            arcpy.Delete_management("DEM_z7")
            arcpy.Delete_management("DEM_z8")
					
        arcpy.AddMessage("*****  All calculations successfuly finished  *****")

		
        return
				
		
		
		
		
		
		
